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ABSTRACT 

The  general  linear  hypothesis  is  formulated  for  a  multivariate 
stationary  stochastic  process.  The  best  (minimum  variance)  linear 
unbiased  estimates  are  derived  for  the  regression  functions  and  it 
is  shown  that  many  signal  estimation  problems  are  special  cases  of 
the  general  linear  model.  Several  examples  are  presented  illustra¬ 
ting  the  technique  for  particular  multivariate  processes. 


INTRODUCTION 


Many  problems  in  the  area  of  applied  time  series  analysis  can 
be  formulated  and  solved  by  extending  and  generalizing  well  known 
techniques  from  the  classical  theory  of  the  multivariate  linear 
hypothesis.  The  analogies  between  some  of  the  physical  models  ex¬ 
pressed  in  terms  of  signals  propagating  across  arrays  and  the  usual 
analysis  of  variance  models  involving  various  kinds  of  effects  are 
striking ,  particular] y  when  it  can  be  assumed  that  the  multivariate 
stochastic  process  in  question  is  weakly  stationary.  In  the  sta¬ 
tionary  case,  simplicity  in  exposition  as  well  as  economy  in  com¬ 
putation  result  from  performing  the  analysis  in  the  frequency  domain 
using  the  properties  of  the  Fourier  transform. 

As  an  example  of  a  simple  estimation  problem  suppose  that  a 
multivariate  stochastic  process  (Y^  (t)  ,  j  =  l,2,...,n,  ~oo<t<oo}  consists 
of  a  signal  which  is  identical  for  each  j  and  a  weakly  stationary 
multivariate  noise  process.  More  specifically  we  assume  that 
Y  .  (t)  =  s(t)  +  n  .  (t)  j  =  l,  2,  .  . .  ,n  (1) 

^  ^  —  oo<  ■£<  oo 

In  this  model  s (t)  is  regarded  as  a  fixed  unknown  signal  with  n^ (t) 
a  weakly  stationary  zero  mean  multivariate  noise  process  with  a 
correlation  function  given  by 


Rjk(t-f)  =  E„.<t)„k(f)  =  C  eiM  <*-*'>  o.k(M)  |5  (2) 

where  we  assume  that  0^(03)  is  the  cross  spectral  density  matrix 
of  the  noise  process.  The  notation  is  simplified  if  we  use  contin¬ 
uous  parameter  processes  and  the  results  apply  equally  well  to  the 
discrete  parameter  sampled  data  version.  By  a  best  linear  unbiased 

A 

estimate  of  the  signal,  say  s(t),  is  meant  the  unbiased  linearly 


filtered  version  of  the  data  (Es  (t)  =  s(t))  which  has  the  smallest 
variance  (E(s(t)  -  s(t))z  =  min.).  This  problem  is  analogous  to 
estimating  the  mean  of  a  set  of  correlated  random  variables  and 
has  been  previously  considered  by  Kelly  and  Levin  5,  and  Capon, 
Greenfield,  and  Kolker,  1.  The  inherent  advantage  of  this  model 
over  the  Wiener  approach  (for  example  see  Papoulis,  8)  is  that  it 
is  not  necessary  to  know  in  advance  the  functional  form  of  the  signal 
in  order  to  minimize  the  mean  square  error  or  variance  of  the  signal 
estimate. 

Frequently  in  the  analysis  of  the  multivariate  process  Y^ (t)  it 
is  necessary  to  consider  more  general  models  for  the  signal  such  as 

(t)  =  s^t)  +  s2  (t-Tj  )  +  n^  (t)  (3) 

where  the  output  of  the  jth  process  is  represented  as  the  sum  of 
two  signals  one  of  which  has  been  delayed  by  a  known  amount  T ^ . 
Equation  (3)  resembles  the  usual  analysis  of  variance  model  with 
Sjl  (t)  corresponding  to  a  row  effect  and  s2(t)  corresponding  to  a 
column  effect.  Delaying  s2  (t)  by  T..  makes  the  separation  possible 
and  corresponds  to  multiplying  by  e  j  in  the  frequency  domain. 
Equation  (3)  suggests  the  utility  of  a  more  general  linear  model 
of  the  form 

P 

yj  (t)  *  C  xjm  Bm(u)  du  +  nj(t)  (4) 

j  =  l, . . .  ,n  — oo<t<«> 

which  includes  p  fixed  unknown  functions  (signals)  (81(t),  P>2(t),..., 
8p  ( t ) )  to  be  estimated  with  Xjk(t)  an  nxp  matrix  of  fixed  functions 
which  are  determined  by  the  model  in  (4).  For  example,  equation 
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(4)  reduces  to  equation  (1)  by  the  choices  ^  (t) 


s(t),  Xjl(t)  * 


6(t)  the  Dirac  delta  function  and  p  =  1.  Equation  (2)  may  be  ob¬ 


tained  by  taking  X.  (t)  =  6ft),  X.  (t)  =  6( t-T . )  with  p  «  2,  and 
(t)  =  (t)  •  $2^  “  s2  *  T^e  Fourier  transforms  of  the  ob¬ 


servables  Xj^(t)  are  often  well  known  for  estimation  problems  and 
the  analysis  proceeds  easily  in  the  frequency  domain.  It  should  be 


noted  that  linear  models  of  the  form  (4)  have  been  developed  for 
the  single  dimensional  case,  Grenander  and  Rosenblatt,  4,  and  the 
multivariate  case  Rosenblatt,  10,  for  some  trigonometric  and  poly¬ 
nomial  regression  models  or  with  special  assumptions  about  the  auto¬ 


regressive  nature  of  the  residuals  n..  (t)  .  Parzen,  9,  also  gives 
the  minimum  variance  solution  for  a  single  dimensional  stochastic 


process.  The  approach  here  takes  advantage  of  formulas  relating 
convolution  in  time  to  multiplication  in  frequency  to  develop  sets 
filters  which  can  be  applied  either  in  time  or  frequency.  In 
addition  we  show  how  the  general  linear  model  can  be  specialized 
to  include  a  majority  of  the  signal  estimation  problems  encounter¬ 
ed  in  practical  situations.  Two  examples  are  presented  using  the 
minimum  variance  unbiased  solutions  for  equations  (3)  and  (4)  when 
it  can  be  assumed  that  the  noises  are  uncorrelated  for  j  /  k. 


We  consider  first  the  derivation  of  the  best  linear  unbiased 
(BLUE)  estimates  for  the  p  regression  functions  (^  (t) 

ep(t))  in  equation  (4)  where  (Xjk(t),  j=l . n  k  =  l,..,,p)  is 

an  nxp  matrix  of  fixed  known  functions.  The  Fourier  transforms  of 
Xjk(t)  and  Bj (t)  are  assumed  to  exist  for  all  j  and  k.  Consider 


(5) 


a  class  of  linear  estimates  of  che  form 

6,<t>  -  l  C  hjk'T>  Yk(t‘T)  dT 

3  k=l  3 

where  h.  (t)  is  a  pxn  matrix  of  filter  functions  to  be  determined. 
3k 

Note  that  determining  h..k(t)  is  sufficient  for  estimating  Bj  (t)  as 
the  filters  need  only  be  convolved  with  the  multivariate  process 

A 

Yj (t)  and  summed  to  produce  Bj (t) . 

To  introduce  the  unbiased  property  we  assume  that 

8  •  (t)  =  [I//h.k(T)Xkin(t-T-u)  Bm  (u)  du  dT  (6) 

3  km-®  J 

which  implies  by  equations  (4)  and  (5)  that 


8  (t)  =  6h  (t)  +  l  C  hik(T)  nk(t_T)  dT 
3  J  k=l  J 


A 

Since  Enk(t)  *  0,  equation  (7)  implies  that  EBj  (t)  =  Bj  (t)  if 
equation  (6)  is  satisfied  for  all  j  and  t.  Then,  by  taking  Fourier 
■transforms  in  equation  (6)  we  arrive  at  the  unbiased  condition  in 
frequency 

ej(t)  -  (8) 

J  mk  J 

where  the  Fourier  transforms  of  the  elements  of  (7)  appear  in  the 
same  order  in  (8).  Now  equation  (8)  implies  that  we  must  have 

J,  ><"»  *  V  <9) 

k=l 

where  <S .  is  the  ordinary  delta  function.  Equation  (9)  can  be 
Dm 

written  in  matrix  notation  as 

HX(oj)  =  I  <10) 

with  I  the  pxp  identity  matrix.  Tne  equation  is  an  identity  at  each 

frequency  w. 

Since  (7)  is  a  vector  of  functions  we  determine  the  estimate 


which  minimizes  the  variance  of 


z(t)  =  l  f  a .  (u)  B .  (t-u)  du 
i=l"  3  3 


an  arbitrary  linearly  filtered  combination  of  the  estimates  6 .  (t) „ 

3 

Note  first  that 

oo 

(z(t)  -  Ez(t))  =  Iff  a .  (u)  h  (t  )  n.  (t-u-r)  du  dT 

j  , k  -co  3  I*  * 

so  that  by  squaring  and  taking  expectations  we  obtain  the  expression 
for  the  variance.  The  final  result  after  taking  Fourier  transforms 
is  expressed  in  the  frequency  domain  as 

E  (z  (t)  —  Ez  (t)  )  2  =  l  l  (  fl  A.  <w)Hjk(aj)okk'  (u)H  (0))  A  *  (<d)  |£ 

j  ,  k  j  1  ,  k  1 


=  /  ah>:h*  a*  (to) 


where  by  A  we  mean  a  lxp  row  vector  consisting  of  the  Fourier  trans¬ 
forms  of  a^ (t)  and  by  A*  we  mean  a  pxl  column  vector  which  is  the 
complex  conjugate  transpose  of  A  with  *  signifying  complex  conjugate. 
Also,  from  equation  (2),  \  =  {Ojk(w)i  is  the  nxn  matrix  of  spectra 
and  cross  spectra  at  the  frequency  u).  In  the  Appendix  we  show  that 
the  inequality 

aB~  1  a*  >  aX  (X*BX)  ~  1X*a*  (13) 

holds  with  B  an  nxn  positive  definite  matrix  and  X  a  pxn  matrix, 
p<n,  of  rank  p  and  a  any  Ixn  complex  row  vector.  Now  with  a  --  AH 
and  \  =  B  1  we  may  derive  a  lower  bound  for  the  variance  given  in 
equation  (12) .  For 


var  z  ( t )  >_  /“  AHX(X*J-1X)_1X*H*A*  (w) 

=  LZ  AfX*!-1^)-1^*  (w) 
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when  the  unbiased  condition  (10)  is  taken  into  account.  Comparing 
(12)  and  (14)  shows  that  equality  is  achieved  when 

H  (cu)  «  (X*^"lX)  ~  1  (to)  (15) 

Hence,  the  time  version  of  the  matrix  of  filters  is  given  by  calcu¬ 
lating  equation  (15)  separately  at  each  ^frequency  and  then  trans¬ 
forming  the  result  back  into  time.  The  resemblance  of  (15)  to  the 
matrix  equations  for  the  usual  weighted  multivariate  hypothesis  is 
obvious . 

We  may  compute  the  variance-covariance  matrix  of  the  estimators 

A  A 

6jL(t)  and  (t)  using  (7).  Again,  squaring,  taking  expectations  and 
transforming  the  results  into  the  frequency  domain  yields 
foov(5.(t),  Sj  CtTj -  c  hZh*  (ui  fs 

-  c  «**-•*>-■  <“>  i  ,i6, 

which  corresponds  to  the  usual  formulas  in  the  multivariate  linear 
hypothesis. 

The  general  solution  given  in  (15)  requires  that  the  spectral 
matrix  of  the  noises  be  known  and  although  this  is  seldom  the  case, 
good  estimates  for  l  can  sometimes  be  obtained  by  using  a  sample  of 

noise  when  the  regression  coefficients  3-(t)  are  not  present.  In 

3 

the  next  section  we  show  how  the  general  linear  model  through  an 
appropriate  choice  for  Xj^(t)  can  be  specialized  to  a  number  of 
interesting  cases  occurring  in  practical  applications. 

SOME  SPECIAL  CASES 

Case  (i)  Estimation  of  the  signal  for  uncorrelated  processes 
Consider  the  model 

Y.  (t)  =  s(t-T.)  +  n.  (t)  (17) 
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which  assumes  that  the  signal  appears  the  same  of  each  process  ex¬ 


cept  that  it  is  delayed  by  a  known  amount  T ^ .  We  assume  also  that 


oo  . 

Rjk  (t- 1 1 )  =  fm  5jko2(w)e 


io)(t-t')  dw 
2? 


which  means  that  the  noise  is  uncorrelated  for  j  /  k  and  has  the 
same  spectrum  (autocorrelation)  at  each  level.  In  this  case  X.^(t) 

»  /i)T 

6  (t-T.. )  and  =  e  j  with  B-^(t)  =  s  (t)  .  Hence  X*  (w)  = 

(e'*'wTl, . .  . ,  e^^n)  and  J  1  (w)  =  (1/o2(uj))  I.  Thus,  equation  (15) 


yields 


H(oj)  =  (l/n)  (eiwri,  eiwT2 , 


so  that 


hk(t)  =  (l/n)  6(t+Tk) 


or  from  (5) 


(t)  =  (l/n)  l  Y  (t+T  ) 


Hence,  the  BLUE  estimate  for  this  case  is  just  the  average  of  the 
process  accounting  for  the  delay  time  Tk  and  the  variance;  of  the 
signal  estimate  is,  from  (16) 

var  s  ( t )  =  /“  o  2  (w)  /n  ^  =  R(0)/n  (22) 

Case  (ii)  Estimation  of  the  signal  for  correlated  processes 

Assume  that  the  basic  model  is 

Yj  (t)  =  s  (t)  +  n^  (t)  (23) 

where  the  correlation  structure  and  power  spectral  matrix  of  the 
noise  is  specified  by  (2).  In  this  case  X^tt)  =  6  (t)  so  that 
Xj^w)  =  1  and  denoting  the  inverse  of  the  spectral  matrix  by  l~ 1 
{ajk(o))}  we  have  by  (15) 
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(24) 


Hk(u) 


Z«jkM 

jfk 


with 


T“;<t»  ’  £  &  (25) 
This  corresponds  to  the  case  considered  in  Kelly  and  Levi„,  5,  and 
Capon,  Greenfield,  and  Kolker,  1. 

Case  (iii)  Estimation  of  multiple  qin„aie  e 

muirxpie  signals  for  uncorrelated 

processes. 

If  in  the  general  model  given  by  equation  (4)  it  is  assumed 
that  the  noise  processes  are  uncorrelated  ^or  J  ,  k  as  in  case  ,i, 
then  IM  I  where  I  is  the  nxn  identity  matrix  and  the 

equations  resemble  those  arising  in  the  multivariate  linear  hypo- 

thesis.  For  (15)  becomes  .  ’ 

» 

Hid))  ='  (X*X)  ~  JX*  (cu)  » 

(26) 

and  (16)  yields 

cov  (8^(t)  (t) )  =  o2(w)  (x*x)_1  (o))^ 

27r  (27) 

AS  a  particular  example  the  multiple  signal  nodel  represented 
by  equation  (3)  yields 


X*  (oj)  = 


i/1/  ...,  1 

eiwTl  _iwT 

1  •  •  •  •  >  e  n 


,x*xr'  <“>  -  mr 


-z  (w) 


-z*(w)  n 
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with 

z(u>)  =  l  e  *wTk  A  (w)  =  n2  -  |  z  (w)  | 2 
k=l 

when  \ (w)  ^  0.  When  A  (w)  =  0  as,  for  example,,  at  u>  =  0,  adjustments 
are  needed  in  the  procedure.  The  next  sections  illustrate  several 
examples  constructed  to  test  the  computational  procedures  and  theory. 
EXAMPLES 

The  theoretical  results  of  the  previous  sections  are  easily 
extended  to  discrete  time  parameter  processes  by  replacing  integrals 
with  summations  and  the  infinite  frequency  range  by  (-■&<&< tt)  but 
in  the  application  to  finite  time  sampled  data  certain  approximations 
are  necessary.  We  must  decide  fi-st  whether  the  analysis  proceeds 
more  reasonably  and  economically  in  the  time  domain  or  the  frequency 
domain.  We  have  chosen  the  frequency  domain  for  several  reasons. 
First,  the  restrictions  imposed  by  signal  models  such  as  equation 
(3)  must  allow  for  time  delays  as  long  as  300  digital  points.  This 
means  that  the  time  domain  analogue  of  the  matrix  product  (X*£-1X)-1 
becomes  extremely  large  if  a  512  point  time  delay  preserving  filter 
is  required.  Hence,  the  frequency  domain  approach  becomes  much 
faster  as  the  analysis  can  be  performed  separately  at  each  frequency. 
In  addition,  for  the  signal  model  discussed  as  Case  (ii) ,  it  has 
been  noted  in  Capon,  Greenfield,  and  Kolker,  1,  that  short  filters 
designed  in  the  time  domain  are  more  sensitive  to  slight  station- 
arity  and  signal  model  perturbations.  Reference  1  also  contains 
several  examples  illustrating  time  and  frequency  domain  computations 
and  concludes  that  the  loss  due  to  using  a  two  sided  infinite  lag 
filter  rather  than  the  physically  realizable  filter  obtained  in  the 
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time  domain  is  small. 

The  approximations  used  in  the  examples  below  are  based  on  the 
use  of  the  finite  Fourier  transform  and  involve  replacing  the  data 
and  filter  by  aliased  versions  of  the  optimum  infinite  two  sided 
Fourier  transforms.  References  2  and  7  contain  excellent  material 
describing  the  theory  and  application  of  the  finite  Fourier  trans¬ 
form.  The  effects  of  aliasing  are  reduced  by  choosing  an  appropriate 
sampling  interval  in  time  based  on  the  highest  frequency  observable 
in  the  data.  In  the  two  examples  given  below  the  signals  are  limited 
to  the  band  from  0  to  10  cycles  per  second  (cps)  with  the  major 
frequency  content  less  than  5  cps.  Hence,  the  data  can  be  sampled 
with  At  =  .05  seconds  yielding  an  aliasing  frequency  of  10  cps. 

The  data  samples  are  long  with  the  first  being  150  seconds  (3,000  pts) 
and  the  second  60  seconds  (1200  points)  so  that  -«><t«»  is  approx¬ 
imately  valid.  Finally,  the  accuracy  of  the  various  approximations 
employed  can  be  evaluated  from  simulated  examples  similar  to  those 
given  below. 

Example  1  (Case  (iii)) 

Here  we  consider  the  model  given  by  equation  (3) 

Yj  (t)  =  s1(t)  +  s^t-Tj)  +  Hj  (t) 

(30) 

j  -  1,2, ... ,20 

Figure  1  shows  the  first  five  channels  of  data  constructed  to 
conform  with  Equation  (30).  Two  seismic  signals  and  smoothed  white 
noise  uncorrelated  from  channel  to  channel  were  added  in  the  same 
proportions  to  produce  twenty  channels.  Five  of  the  twenty  channels 
are  displayed  in  Figure  1.  The  general  signal  estimation  procedure 
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was  to  compute  the  matrix  product  (X*X)_,X*  (u>)  at  each  of  256  fre¬ 
quencies  using  Equations  (28)  and  (29).  The  singularity  at  w-0  was 
eliminated  by  tapering  the  frequency  response  functions  down  to  zero 
at  0  cps.  A  fast  Fourier  transform  subroutine  (see  McCowan,7)  applied 
to  the  frequency  response  matrix  H(to)  produced  the  impulse  response 
functions  shown  in  Figure  2.  Note  that  the  net  results  of  the  filters 
are  to  reinforce  and  sum  the  aligned  signal  s^t)  while  canceling  the 
second  signal  s2(t)  which  appears  at  the  given  time  delays.  The  re¬ 
verse  holds  when  estimating  s2(t)  with  the  unaligned  signal  reinforced 
and  the  aligned  signal  canceled.  The  signal  estimates  are  formed  using 
a  high  speed  convolution  subroutine  which  calculates  equation  (5)  as 
a  matrix  product  in  the  frequency  domain  and  then  re-transforms  the 
result  to  obtain  time  domain  estimates  for  the  regression  functions. 
Note  that  in  the  application  of  the  convolution  each  Y^  (t)  is  avail¬ 
able  only  for  the  finite  interval  T  and  We  assume  that  Yj (t)  =  0  for 
t  X  T.  The  true  signals  s^t)  and  s2  (t)  are  displayed  along  with 
the  estimates  s1  (t)  and  s2(t)  in  Figure  1.  The  generation  of  the 
filters  at  256  frequencies  or  512  time  points  for  20  channels  took 
eight  minutes  of  CDC  1604  time  with  the  subsequent  convolution  re¬ 
quiring  thirty  minutes.  However,  by  combining  the  operations  in 

a  different  order  the  total  running  time  has  been  reduced  to  eight 
minutes. 


Example  2 


As  an  example  of  a  more  complicated  model  consider  the  multi¬ 
variate  process  shown  in  Figura  3.  Each  level  or  channel  conforms 
artificially  to  the  model. 
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(30) 


Y,(t)  -  B,(t)  +  1  C  X.  (t-u)  3  (u)  du  +  n .  (t) 

J  A  m*2  J  J 


(j  *  1#  .  .  .  /  5) 

The  data  was  generated  by  adding  noise  and  the  known  functions  3j  (t) 
in  Figure  3  to  the  data  using  the  transfer  functions  Xj^(t).  The 
real  valued  frequency  responses  (X^U),  f  =  w/2ir  cps)  are  shown  in 
Figure  4  for  the  range  0  to  2  cps.  The  response  functions  were  low 
pass  filtered  with  a  cutoff  at  2  cps:  hence,  there  is  no  signi¬ 
ficant  frequency  response  above  that  value.  The  optimum  filters  in 
this  case  were  200  points  or  10  seconds  long  and  were  computed  by 
evaluating  the  matrix  product  (X*X)  1 X*  )  at  200  frequencies.  (The 
fast  Fourier  transform  was  not  used  here  so  the  number  of  data  points 
is  not  a  power  of  two.)  The  resulting  matrix  of  filter  coefficients 
was  convolved  with  the  multivariate  process  using  equation  (5)  to 
generate  the  estimates  shown  in  Figure  3.  Comparison  of  these  es¬ 
timated  regression  functions  with  the  true  regression  functions  in¬ 
dicates  that  the  procedure  again  produces  reasonable  estimates. 
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APPENDIX 


PROOF  OF  BASIC  INEQUALITY 

In  order  to  prove  the  inequality  (13)  we  use  the  n  variate 
complex  normal  distribution  defined  by  Goodman,  3.  An  n-variate 
complex  normal  random  variable  =  (?^»  ?2'  an  n_tuP^e 

of  complex  random  variables  with  a  probability  density  function 

given  by 

(£)  =  |tt£|  n  exp  -  (^  -  a)B  1  (£  -  a )  *  (31) 

with 

B  =  E(C  -  a)  *  (C  -  a) 

and  B  a  Hermitian  positive  definite  complex  covariance  matrix.  The 
mean  of  the  complex  vector  K  is  assumed  to  be  the  complex  vector  a. 

We  shall  need  to  calculate  the  discrimination  information 
measure  as  developed  in  Fullback,  6.  If  by  -*-s  mean^  a  zero 

mean  complex  multivariate  normal  random  variable  with  covariance 
matrix  B,  the  discrimination  information  is  defined  to  be 

p 

I  (1: 2; 5)  =  I  P,  (S)  log  — 1  ^  d?  (32) 

-  5  1  ~ 

which  for  Pj^  (O  and  P2(0  above  becomes 

I (Is  2; O  =  a  B-1  a*  (33) 

Now  if  X  is  an  nxp  (p<n)  matrix  of  rank  p;  define  the  transformation 
n  =  £X  (34) 

with  the  pxl  vector  n_  having  a  complex  multivariate  normal  distri¬ 
bution  with  mean  aX  and  variance-covariance  matrix  X*BX.  Then  X*BX 
is  non-singular  and  Hermitian  positive  definite  and 

I  (1:2;  n)  =  aX  (X*BX)  ”  1  X*a*  (35) 
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I 


Now  since  (34)  is  a  measurable  transformation  we  have  (for  the  proof 
of  this  inequality  in  the  real  ca'ie,  see  Kullback,  6,  p.57)  I(li2;r[) 
£  I(l:2;£)  which  implies  the  inequality 

aB-  1_o*  >  aX(X*BX)_1X*a*  (36) 
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Figure  1.  Separation  of  Two  Mixed  Signals 
By  Minimum  Variance  Filters 
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if  Observations  X..  (f),  cj=2Trf 


